Pre-processing

counts <- read.csv('../data/geo_submission/counts_raw.csv', sep = ',', header = T,row.names = 1)
coldata <- read.csv('../data/geo_submission/cell_metadata_raw.csv', row.names = 1) 

# identify spike-in genes
ercc_rows <- grepl('^ERCC', rownames(counts))
ercc_counts <- counts[ercc_rows,]

# remove spike-in genes from counts
counts <- counts[!ercc_rows,]

cell filtering

Plate P4rep exhibits a large number of poor-quality cells, known to invalidate the outlier-detection approach for cell-filtering. To obviate, QC metrics are shared across plates to select high-quality cells.

mito <- grep('^mt.*', rowData(sce)[,'SYMBOL'])

stats <- perCellQCMetrics(sce, subsets=list(Mt=mito))

# share QC metrics across good plates
discard_libsize <- isOutlier(stats$sum, type = 'lower', log = T, 
                             batch = sce$plate, subset = sce$plate %in% c('2','7'))
discard_features <- isOutlier(stats$detected, type = 'lower', log = T, 
                              batch = sce$plate, subset = sce$plate %in% c('2','7'))
discard_mito <- isOutlier(stats$subsets_Mt_percent, type = 'higher', log = F, 
                          batch = sce$plate, subset = sce$plate %in% c('2','7'))
discard_ercc <- isOutlier(stats$altexps_spike_in_percent, type = 'higher', log = F,
                          batch = sce$plate, subset = sce$plate %in% c('2','7'))
discard_df <- data.frame(libsize = discard_libsize, features = discard_features, 
                         mitoc = discard_mito, ercc = discard_ercc, plate = sce$plate)

# discard cells with low library size / few detected genes / high ercc counts
discard_df$discard <- sign(rowSums(discard_df[,c(1,2,4)]))

colSums(discard_df[,c('libsize','features','mitoc','ercc','discard')])
##  libsize features    mitoc     ercc  discard 
##      182      228      109      253      322
colData(sce)[,colnames(stats)] <- stats
sce$discard <- as.logical(discard_df$discard)

qc_plot <- wrap_plots(
    plotColData(sce, x="plate", y="sum", colour_by="discard", point_size =.25,
                show_median =T) +
      scale_y_log10() + labs(y = 'library size'),
    plotColData(sce, x="plate", y="detected", colour_by="discard", point_size =.25,
                show_median =T) + scale_y_log10() + labs(y = 'detected genes'),
    plotColData(sce, x="plate", y="subsets_Mt_percent", colour_by="discard",
                point_size =.25, show_median =T) +  labs(y = 'mitochondrial_pct'),
    plotColData(sce, x="plate", y="altexps_spike_in_percent", colour_by="discard",
                point_size =.25, show_median =T) + ggtitle("spike_in percent"))
qc_plot <- qc_plot  + plot_layout(guides = 'collect') & theme_bw(base_size = 7)

print(qc_plot)

# discard low-quality cells
sce <- sce[,!sce$discard]

Normalization via deconvolution Pool-based normalization is the recommended normalization method to remove bias associated with RNA content [Lun, A. et al., 2017]. Here, a version optimized for handling batch effects is used:

sce <- computeSumFactors(sce)
sce <- multiBatchNorm(sce, batch = sce$plate)

Feature selection The mean-variance relationship is modelled separately for each plate, and significant genes are retained:

# use spike-in to model gene-variance relationship
genevar <- modelGeneVarWithSpikes(sce, "spike_in", block=sce$plate)
chosen.hvgs <- getTopHVGs(genevar, fdr.threshold = 0.05)
# remove cell-cycle genes 
cc.genes <- AnnotationDbi::select(org.Mm.eg.db, keys="GO:0007049", keytype="GOALL", column="ENSEMBL")$ENSEMBL
chosen.hvgs <- chosen.hvgs[!rowData(sce)[chosen.hvgs, 'ENSEMBL'] %in% cc.genes]
genevar_plots = lapply(1:length(unique(sce$plate)), function(i){
  genevar_genes =  genevar$per.block[[i]] %>% as.data.frame()
  genevar_spikeins = data.frame(mean = metadata(genevar$per.block[[i]])$mean, 
                                var = metadata(genevar$per.block[[i]])$var)
  genevar_genes$trend = metadata(genevar$per.block[[i]])$trend(genevar_genes$mean)
  ggplot(genevar_genes, aes(x = mean, y = total)) + geom_point(size = 0.25*pointsize, col = 'grey60') +
    geom_point(inherit.aes = F, data = genevar_spikeins, aes(x = mean, y = var),size = pointsize, col = 'red') +
    geom_line(inherit.aes = F, data = genevar_genes, aes(x = mean, y = trend), col = 'blue') +
    theme_bw() + labs(title = unique(sce$plate)[i]) + xlim(c(0,12)) + ylim(c(0,17.5))
})

wrap_plots(genevar_plots)

batch correction Plate-effects can be corrected by using a linear method from the limma package:

assay(sce, "corrected") <- removeBatchEffect(logcounts(sce), batch = sce$plate)

Analysis

dimensionality reduction

sce <- runPCA(sce, exprs_values = "corrected", subset_row = chosen.hvgs)
pc_var = data.frame(var_pct = attr(reducedDim(sce),'percentVar'))
ggplot(data.frame(var_pct = attr(reducedDim(sce),'percentVar')), aes(x = order(var_pct, decreasing = T), y = var_pct)) + 
  geom_point()+ geom_line() + theme_bw() + geom_vline(xintercept = 10, col = 'red')

Analysis of explained variance per PC shows that most information is contained in the first few components. We will retain 10 for downstream analyses.

PHATE has been popularized as an effective dimensionality reduction method to visualize differentiation trajectories:

reducedDim(sce, 'PHATE2D') <-  phateR::phate(reducedDim(sce,'PCA')[,1:10], ndim =2, verbose =0)$embedding
coldata_dr(sce,color.column = 'cell_type2', dimred = 'PHATE2D',size = pointsize)  +
  scale_color_tableau(palette = 'Tableau 10') + labs( title = 'cell type')

This visualization broadly confirms proximity relationships across hematopoietic compartments.

diffusion pseudotime and lineage scores

Additional annotation is generated by:

  • Computing pseudotime coordinates using diffusion pseudotime via the destiny package.

  • Running Seurat’s AddModuleScore function to aggregate the expression of known lineage markers and score their expression on each cell:

# compute diffusion map (necessary for diffusion pseudotime)
dmap = destiny::DiffusionMap(reducedDim(sce,'PCA')[,1:10])
reducedDim(sce, 'DiffusionMap') <-  dmap@eigenvectors

# isolate stem cell tip on PHATE coordinates
stem_tip = as.integer(which.min(as.logical(sce$cell_type2 == 'HSC') * reducedDim(sce,'PHATE2D')[,1]))
sce$stem_tip = FALSE
sce$stem_tip[stem_tip] = TRUE

dpt =destiny::DPT(dmap, tips = stem_tip)

# assign pseudotime coordinate
sce$dpt = dpt$DPT489 # this is the value of stem_tip
sce$dpt_rank = rank(sce$dpt)

lineage_sets <- openxlsx::read.xlsx('../docs/camargo_sets.xlsx')

# load SingleCellExperiment object in Seurat
seu <- CreateSeuratObject(counts = counts(sce))
seu@assays$RNA@data = assay(sce, 'corrected')
selection = c('preMeg','preEry','preGM','preB')

# compute scores for selected sets
for(i in selection){
  feature_names = rownames(seu)[rowData(sce)[, 'SYMBOL'] %in% lineage_sets[[i]]]
  seu <- AddModuleScore(seu, features = list(feature_names), name = i)}

# transfer scores in SingleCellExperiment object
colData(sce)[,selection] <- seu@meta.data[,paste0(selection,'1')]

# plotting
p_dpt <- coldata_dr(sce, color.column = 'dpt_rank', dimred = 'PHATE2D',size = pointsize) +
  scale_color_viridis_c() + guides(col = F) + labs( title ='pseudotie rank')
p_my <- coldata_dr(sce, color.column = 'preGM', dimred = 'PHATE2D',size = pointsize) + 
  scale_color_viridis_c()  + guides(col = F) + labs( title ='myeloid score')
p_mk <- coldata_dr(sce, color.column = 'preMeg', dimred = 'PHATE2D',size = pointsize) +
  scale_color_viridis_c()  + guides(col = F) + labs( title ='mk score')
p_ery <- coldata_dr(sce, color.column = 'preEry', dimred = 'PHATE2D',size = pointsize) +
  scale_color_viridis_c()  + guides(col = F) + labs( title ='ery score')

fig <- (p_dpt +   p_mk + p_my + p_ery) + plot_layout(guides = 'collect') & 
  theme(axis.title = element_blank(),
        axis.text = element_blank())

fig

Next, we’ll investigate proximity relations between HSCs and MkP subsets.

proximity analysis

MkP_CD48 vs Sca1

Distance between CD48 subsets of MkP and HSCs:

mkp_cols <- sce$cell_type2 == 'MkP'
hsc_cols <- sce$cell_type2 == 'HSC' 
distance <- as.matrix(dist(reducedDim(sce)[,1:10]))
distance <- distance[mkp_cols, hsc_cols] %>% rowMedians()

df <- data.frame(subset = sce$CD48_binary[mkp_cols], distance = distance)

distance_boxplot_MkP <- ggboxplot(data = df, x = 'subset', y = 'distance', add = 'jitter', color = 'subset', palette = viridis(2)) + 
  stat_compare_means(comparisons = list(c('low','high')), 
                     method = 'wilcox.test',label = "p.signif") + xlab('MkP CD 48 subset') + ylab('proximity to HSC') + scale_y_reverse() + 
  guides(col =F)

distance_boxplot_MkP

HSC_Sca1 subsets vs MkP_48low

Distance between Sca-1 subsets of HSCs and CD48-/lo MkPs:

mkp_cols <- sce$cell_type2 == 'MkP' & sce$CD48_binary == 'low'
hsc_cols <- sce$cell_type2 == 'HSC' 
distance <- as.matrix(dist(reducedDim(sce)[,1:10]))
distance <- distance[mkp_cols, hsc_cols] %>% colMedians()

df <- data.frame(subset = sce$Sca_1_binary[hsc_cols], distance = distance, Sca1 = sce$Sca_1_log10[hsc_cols], 
                 CD201 = sce$CD201_log10_normalized[hsc_cols])

distance_boxplot_HSC <- ggboxplot(data = df, x = 'subset', y = 'distance', add = 'jitter', color = 'subset', palette = viridis(2)) + 
  stat_compare_means(comparisons = list(c('low','high')), 
                     method = 'wilcox.test',label = "p.signif") + xlab('HSC Sca1 subset') + ylab('proximity to MkP (CD 48 low)') + scale_y_reverse() + 
  guides(col =F)

distance_boxplot_HSC

##

refine populations

Based on expression of surface markers, we refine the populations. Subsets that we’re interested in are higlighted in the following plot:

Note: Very few cells are observed as Sca1- and CD201+. Here, CD201+ HSCs are grouped for simplicity.

sce$cell_type3<- as.character(sce$cell_type2)

sce$cell_type3[sce$cell_type2 == 'MkP' & sce$CD48_binary == 'low'] = 'MkP_CD48-'
sce$cell_type3[sce$cell_type2 == 'MkP' & sce$CD48_binary == 'high'] = 'MkP_CD48+'

sce$cell_type3[sce$cell_type2 == 'HSC' & sce$Sca_1_binary == 'low'] = 'HSC_Sca1-'
sce$cell_type3[sce$cell_type2 == 'HSC' & sce$Sca_1_binary == 'high'] = 'HSC_Sca1+'

sce$cell_type3[sce$cell_type2 == 'MPP' & sce$Sca_1_binary == 'low'] = 'MPP_Sca1-'
sce$cell_type3[sce$cell_type2 == 'MPP' & sce$Sca_1_binary == 'high'] = 'MPP_Sca1+'

sce$cell_type3[sce$cell_type2 == 'HPC_1' & sce$Sca_1_binary == 'low'] = 'HPC1_Sca1-'
sce$cell_type3[sce$cell_type2 == 'HPC_1' & sce$Sca_1_binary == 'high'] = 'HPC1_Sca1+'

sce$cell_type3[sce$cell_type3 == 'HSC_Sca1+' & sce$CD201_binary == 'high'] = 'HSC_Sca1+__CD201+'
sce$cell_type3[sce$cell_type3 == 'HSC_Sca1+' & sce$CD201_binary == 'low'] = 'HSC_Sca1+__CD201-'

sce$cell_type3[sce$cell_type3 == 'HSC_Sca1-' & sce$CD201_binary == 'high'] = 'HSC_Sca1-__CD201+'
sce$cell_type3[sce$cell_type3 == 'HSC_Sca1-' & sce$CD201_binary == 'low'] = 'HSC_Sca1-__CD201-'

# transcript-based
sce$cell_type3[sce$cell_type2 == 'HPC_1' & sce$preB >2] = 'HPC1_lymphoid'
sce$cell_type4<- as.character(sce$cell_type3)

sce$cell_type4[sce$cell_type3 == 'HSC_Sca1-__CD201+'] = 'HSC_CD201+'
sce$cell_type4[sce$cell_type3 == 'HSC_Sca1+__CD201+'] = 'HSC_CD201+'
p <- coldata_dr(sce, col = 'cell_type4',dimred = 'PHATE2D', size = pointsize) +
  scale_color_manual(values = c(rep('lightgrey',3),"#4E79A7","#F28E2B","#E15759", rep('lightgrey',2),"#00AFBB", "#E7B800",rep('lightgrey',2)))  + labs(title = 'hsc subsets')

p

PAGA analysis

Proximity analysis for n>2 populations is unwieldy. We will use PAGA to investigate proximity relations systematically:

### assign gene expression values in easily accessible containers:
pca <- reducedDim(sce,'PCA')
logc <- assay(sce,'corrected')
coldata <- as.data.frame(colData(sce))
rowdata <- as.data.frame(rowData(sce))
import scanpy as sc
import numpy as np
x = sc.AnnData(X = r.logc.transpose())
x.obsm['X_pca'] = np.array(r.pca)
x.obs = r.coldata
x.var_names = r.rowdata.index

sc.pp.neighbors(x, n_pcs =10)
sc.tl.paga(x, groups='cell_type4')
## /Users/greco/opt/anaconda3/lib/python3.7/site-packages/anndata/_core/anndata.py:1192: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
##   if is_string_dtype(df[key]) and not is_categorical(df[key])
## ... storing 'plate' as categorical
## ... storing 'mouse' as categorical
## ... storing 'age' as categorical
## ... storing 'sex' as categorical
## ... storing 'cell_type' as categorical
## ... storing 'LR' as categorical
## ... storing 'Row' as categorical
## ... storing 'cell_type2' as categorical
## ... storing 'CD48_binary' as categorical
## ... storing 'Sca_1_binary' as categorical
## ... storing 'CD201_binary' as categorical
## ... storing 'cell_type3' as categorical
## ... storing 'cell_type4' as categorical
sc.pl.paga(x, color=['cell_type4'], threshold =.9, layout = 'fr', edge_width_scale =0.25)

The PAGA layout largely confirms the results from the labelling propagation data. To verify that this network is the “best fitting one”, i.e., no alternative models can explain the as well as the one displayed, we can plot the distribution of the confidence placed by the model in each link, and highlight the links that exceed the selected threshold:

connectivies = x.uns['paga']['connectivities']
connectivities = as.matrix(py$connectivies)
connectivities = unlist(connectivities[upper.tri(connectivities)])
df = data.frame(confidence = connectivities)

fig_sx2 = ggplot(df, aes(x = confidence)) + geom_histogram(aes(fill = confidence > 0.9),bins = 20) + 
  scale_fill_tableau() + theme(text = element_text(size =10))
fig_sx2

Pseudotime vs Mk score

PAGA connectivities seem to infer a differentiaion link from CD48- MkP to CD48+ MkP. By looking more closely at the relationship between pseudotime and Mk score, we observe how differentiation (of which Mk score represents a transcriptional proxy) proceeds indipendently in the two subsets:

# use default dpt
df = colData(sce)[,c('cell_type3','cell_type4','preMeg','dpt', 'CD48_log10')] %>%
  as.data.frame() %>%
  dplyr::filter(grepl('^MkP', cell_type3), preMeg > 0)

p2 <- ggscatter(df%>% dplyr::filter(grepl('^MkP', cell_type3)), x = 'dpt', y = 'preMeg', color = "cell_type3",
                size =pointsize,
                add = "reg.line",  # Add regressin line
                #add.params = list(color = 'cell_type3'),
                palette =c("#00AFBB", "#E7B800"),
                add.params = list(color = "cell_type3", fill = "lightgray"), # Customize reg. line
                conf.int = TRUE, # Add confidence interval
                cor.coeff.args = list(label.sep = "\n")) + labs(y = 'mk score', x = 'pseudotime') + 
  stat_cor(aes(color = cell_type3), size =10/.pt) + guides(col='none')

fig_sx3 = p2 + plot_annotation(caption = 'two cells from myeloid branch with low mk score were excluded') 
fig_sx3
## `geom_smooth()` using formula 'y ~ x'

Other plots:

Sca-1, CD48 surface expression

df <- as.data.frame(reducedDim(sce, 'PHATE2D')[,c(1,2)])
colnames(df) <- c('dr1','dr2')
df[['CD48_log10']] <- colData(sce)[,'CD48_log10']
df[['Sca1_log10']] <- colData(sce)[,'Sca_1_log10']

p_Sca1 <- ggplot(df, aes(x = dr1, y = dr2)) + geom_point(col = 'grey90', alpha =.8, size =pointsize,stroke =0) + 
  geom_point(data= df[sce$cell_type2 %in% c('HSC','MPP','HPC_1'),] , aes(x = dr1,y=dr2, col = Sca1_log10), size =pointsize, stroke=0,inherit.aes = F)  + scale_color_viridis_c() + guides(col =F) + labs( title ='Sca1')

p_mkpsubs <- coldata_dr(sce, col = 'cell_type3',dimred = 'PHATE2D', size = pointsize) +
  scale_color_manual(values = c(rep('lightgrey',9),"#4E79A7","#F28E2B", rep('lightgrey',2)))  + labs(title = 'mkp subsets')

fig_3ef <- (p_Sca1 / p_mkpsubs) & 
  theme(legend.position = 'none',
        axis.title = element_blank(),
        axis.text = element_blank())
fig_3ef

proximity analysis in PHATE space

MkP_CD48 vs Sca1

mkp_cols <- sce$cell_type2 == 'MkP'
hsc_cols <- sce$cell_type2 == 'HSC' 
distance <- as.matrix(dist(reducedDim(sce,'PHATE2D')))
distance <- distance[mkp_cols, hsc_cols] %>% rowMedians()

df <- data.frame(subset = sce$CD48_binary[mkp_cols], distance = distance, 
                 Sca1 = sce$Sca_1_log10[mkp_cols], 
                 CD201 = sce$CD201_log10_normalized[mkp_cols])

distance_boxplot_MkP <- ggboxplot(data = df, x = 'subset', y = 'distance', add = 'jitter', 
                                  color = 'subset', add.params = list(size = pointsize), palette = viridis(2)) + 
  stat_compare_means(comparisons = list(c('low','high')), 
                     method = 'wilcox.test',label = "p.signif") + xlab('MkP CD48 subset') + ylab('proximity to HSC') + scale_y_reverse() + guides(col =F)
distance_boxplot_MkP

HSC_Sca1 subsets vs MkP_48low

mkp_cols <- sce$cell_type2 == 'MkP' & sce$CD48_binary == 'low'
hsc_cols <- sce$cell_type2 == 'HSC' 
distance <- as.matrix(dist(reducedDim(sce,'PHATE2D')))
distance <- distance[mkp_cols, hsc_cols] %>% colMedians()

df <- data.frame(subset = sce$Sca_1_binary[hsc_cols], distance = distance, Sca1 = sce$Sca_1_log10[hsc_cols], 
                 CD201 = sce$CD201_log10_normalized[hsc_cols])

distance_boxplot_HSC <- ggboxplot(data = df, x = 'subset', y = 'distance', add = 'jitter', 
                                  color = 'subset', add.params = list(size = pointsize), palette = viridis(2)) + 
  stat_compare_means(comparisons = list(c('low','high')), 
                     method = 'wilcox.test',label = "p.signif") + xlab('HSC Sca1 subset') + ylab('proximity to MkP (CD 48 low)') + scale_y_reverse() + guides(col =F)


fig_sx1 <- (distance_boxplot_MkP / distance_boxplot_HSC) & theme(text = element_text(size =7), axis.text.y = element_blank())

fig_sx1

Vwf and CD41 expression across cell_types

my_theme <- theme(axis.text = element_blank(),
                  axis.title = element_blank()) 
p_vwf <- gene_dr(sce, gene = 'Vwf', dimred = 'PHATE2D', size = pointsize) + guides(col =F) + my_theme
p_cd41 <- gene_dr(sce, gene = 'Itga2b', dimred = 'PHATE2D', size = pointsize) + guides(col =F) + my_theme

# boxplots -------------------------------------------------------

cell_types = c('HSC_CD201+', 'HSC_Sca1+__CD201-', 'HSC_Sca1-__CD201-',  'MkP_CD48-', 'MkP_CD48+','HPC1_Sca1+')

assayname <- 'logcounts'

df <- colData(sce)[sce$cell_type4 %in% cell_types,] %>% as.data.frame() %>%
  dplyr::select(cell_type4, plate,CD41_log10 )
df$Vwf <- assay(sce, assayname)['Vwf', sce$cell_type4 %in% cell_types]
df$Cd41 <- assay(sce, assayname)['Itga2b', sce$cell_type4 %in% cell_types]
df$cell_type4 <- factor(df$cell_type4, 
                        levels = c('HSC_CD201+', 'HSC_Sca1+__CD201-', 'HSC_Sca1-__CD201-',  'MkP_CD48-', 'MkP_CD48+','HPC1_Sca1+'))

comparisons <-  list(c('HSC_CD201+', 'HSC_Sca1+__CD201-'),
                     c('HSC_CD201+', 'HSC_Sca1-__CD201-'),
                     c('HSC_CD201+', 'MkP_CD48-'),
                     c('HSC_CD201+', 'MkP_CD48+'),
                     c('HSC_Sca1+__CD201-', 'HSC_Sca1-__CD201-'),
                     c('HSC_Sca1+__CD201-', 'MkP_CD48-'),
                     c('HSC_Sca1+__CD201-', 'MkP_CD48+'),
                     c('HSC_Sca1-__CD201-', 'MkP_CD48-'),
                     c('HSC_Sca1-__CD201-', 'MkP_CD48+'),
                     c('MkP_CD48-', 'MkP_CD48+'))



p_vwf.2 <- ggplot(df, aes(x = cell_type4, y = Vwf)) + geom_boxplot(outlier.colour = NA) + 
   geom_jitter(stroke =0 , alpha =0.5, size = pointsize) + stat_compare_means(comparisons = comparisons, label = 'p.signif', size =7/.pt) +
  theme(axis.text.x = element_text(angle = -30, hjust =0, vjust = 0.5)) + theme(axis.text.x = element_blank())

p_cd41.2 <- ggplot(df, aes(x = cell_type4, y = Cd41)) + geom_boxplot(outlier.colour = NA) + 
   geom_jitter(stroke =0 , alpha =0.5, size = pointsize) + stat_compare_means(comparisons = comparisons, label = 'p.signif', size =7/.pt) +
  theme(axis.text.x = element_text(angle = -30, hjust =0, vjust = 0.5))

fig_s7e <- (p_vwf + p_vwf.2 + p_cd41 + p_cd41.2) + 
  plot_layout(ncol =2, widths = c(0.4,0.8)) & theme(text = element_text(size =7))
fig_s7e

sessionInfo - R and Python

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ensembldb_2.14.1            AnnotationFilter_1.14.0     GenomicFeatures_1.42.3     
##  [4] viridis_0.6.1               viridisLite_0.4.0           patchwork_1.1.1            
##  [7] ggpubr_0.4.0                ggthemes_4.2.4              forcats_0.5.1              
## [10] stringr_1.4.0               dplyr_1.0.6                 purrr_0.3.4                
## [13] readr_1.4.0                 tidyr_1.1.3                 tibble_3.1.2               
## [16] tidyverse_1.3.1             limma_3.46.0                org.Mm.eg.db_3.12.0        
## [19] AnnotationDbi_1.52.0        AnnotationHub_2.22.1        BiocFileCache_1.14.0       
## [22] dbplyr_2.1.1                batchelor_1.6.3             SeuratObject_4.0.2         
## [25] Seurat_4.0.2                scater_1.18.6               ggplot2_3.3.3              
## [28] scran_1.18.7                SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
## [31] Biobase_2.50.0              GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
## [34] IRanges_2.24.1              S4Vectors_0.28.1            BiocGenerics_0.36.1        
## [37] MatrixGenerics_1.2.1        matrixStats_0.59.0          reticulate_1.20            
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3                rtracklayer_1.50.0           
##   [3] scattermore_0.7               bit64_4.0.5                  
##   [5] knitr_1.33                    irlba_2.3.3                  
##   [7] DelayedArray_0.16.3           data.table_1.14.0            
##   [9] rpart_4.1-15                  RCurl_1.98-1.3               
##  [11] generics_0.1.0                cowplot_1.1.1                
##  [13] RSQLite_2.2.7                 RANN_2.6.1                   
##  [15] proxy_0.4-26                  future_1.21.0                
##  [17] bit_4.0.4                     xml2_1.3.2                   
##  [19] lubridate_1.7.10              spatstat.data_2.1-0          
##  [21] httpuv_1.6.1                  assertthat_0.2.1             
##  [23] xfun_0.23                     hms_1.1.0                    
##  [25] jquerylib_0.1.4               evaluate_0.14                
##  [27] promises_1.2.0.1              progress_1.2.2               
##  [29] DEoptimR_1.0-9                fansi_0.5.0                  
##  [31] readxl_1.3.1                  igraph_1.2.6                 
##  [33] DBI_1.1.1                     htmlwidgets_1.5.3            
##  [35] spatstat.geom_2.1-0           ellipsis_0.3.2               
##  [37] RSpectra_0.16-0               backports_1.2.1              
##  [39] biomaRt_2.46.3                deldir_0.2-10                
##  [41] sparseMatrixStats_1.2.1       vctrs_0.3.8                  
##  [43] TTR_0.24.2                    ROCR_1.0-11                  
##  [45] abind_1.4-5                   cachem_1.0.5                 
##  [47] RcppEigen_0.3.3.9.1           withr_2.4.2                  
##  [49] robustbase_0.93-8             vcd_1.4-8                    
##  [51] sctransform_0.3.2             GenomicAlignments_1.26.0     
##  [53] prettyunits_1.1.1             xts_0.12.1                   
##  [55] goftest_1.2-2                 cluster_2.1.2                
##  [57] lazyeval_0.2.2                laeken_0.5.1                 
##  [59] crayon_1.4.1                  labeling_0.4.2               
##  [61] edgeR_3.32.1                  pkgconfig_2.0.3              
##  [63] ProtGenerics_1.22.0           nlme_3.1-152                 
##  [65] vipor_0.4.5                   nnet_7.3-16                  
##  [67] rlang_0.4.11                  globals_0.14.0               
##  [69] lifecycle_1.0.0               miniUI_0.1.1.1               
##  [71] modelr_0.1.8                  rsvd_1.0.5                   
##  [73] cellranger_1.1.0              polyclip_1.10-0              
##  [75] RcppHNSW_0.3.0                lmtest_0.9-38                
##  [77] Matrix_1.3-4                  carData_3.0-4                
##  [79] boot_1.3-28                   zoo_1.8-9                    
##  [81] reprex_2.0.0                  beeswarm_0.4.0               
##  [83] ggridges_0.5.3                png_0.1-7                    
##  [85] knn.covertree_1.0             bitops_1.0-7                 
##  [87] KernSmooth_2.23-20            Biostrings_2.58.0            
##  [89] blob_1.2.1                    DelayedMatrixStats_1.12.3    
##  [91] parallelly_1.26.0             rstatix_0.7.0                
##  [93] ggsignif_0.6.1                beachmat_2.6.4               
##  [95] scales_1.1.1                  memoise_2.0.0                
##  [97] magrittr_2.0.1                plyr_1.8.6                   
##  [99] hexbin_1.28.2                 ica_1.0-2                    
## [101] zlibbioc_1.36.0               compiler_4.0.4               
## [103] dqrng_0.3.0                   RColorBrewer_1.1-2           
## [105] pcaMethods_1.82.0             fitdistrplus_1.1-5           
## [107] Rsamtools_2.6.0               cli_2.5.0                    
## [109] XVector_0.30.0                listenv_0.8.0                
## [111] pbapply_1.4-3                 ggplot.multistats_1.0.0      
## [113] MASS_7.3-54                   mgcv_1.8-36                  
## [115] tidyselect_1.1.1              phateR_1.0.7                 
## [117] stringi_1.6.2                 highr_0.9                    
## [119] yaml_2.2.1                    askpass_1.1                  
## [121] BiocSingular_1.6.0            locfit_1.5-9.4               
## [123] ggrepel_0.9.1                 grid_4.0.4                   
## [125] sass_0.4.0                    tools_4.0.4                  
## [127] future.apply_1.7.0            rio_0.5.26                   
## [129] rstudioapi_0.13               bluster_1.0.0                
## [131] foreign_0.8-81                gridExtra_2.3                
## [133] smoother_1.1                  farver_2.1.0                 
## [135] scatterplot3d_0.3-41          Rtsne_0.15                   
## [137] digest_0.6.27                 BiocManager_1.30.15          
## [139] shiny_1.6.0                   Rcpp_1.0.6                   
## [141] car_3.0-10                    broom_0.7.6                  
## [143] scuttle_1.0.4                 BiocVersion_3.12.0           
## [145] later_1.2.0                   RcppAnnoy_0.0.18             
## [147] httr_1.4.2                    colorspace_2.0-1             
## [149] XML_3.99-0.6                  rvest_1.0.0                  
## [151] fs_1.5.0                      tensor_1.5                   
## [153] ranger_0.12.1                 splines_4.0.4                
## [155] uwot_0.1.10                   statmod_1.4.36               
## [157] spatstat.utils_2.1-0          sp_1.4-5                     
## [159] plotly_4.9.4                  xtable_1.8-4                 
## [161] jsonlite_1.7.2                destiny_3.4.0                
## [163] R6_2.5.0                      pillar_1.6.1                 
## [165] htmltools_0.5.1.1             mime_0.10                    
## [167] glue_1.4.2                    fastmap_1.1.0                
## [169] VIM_6.1.0                     BiocParallel_1.24.1          
## [171] BiocNeighbors_1.8.2           class_7.3-19                 
## [173] interactiveDisplayBase_1.28.0 codetools_0.2-18             
## [175] utf8_1.2.1                    lattice_0.20-44              
## [177] bslib_0.2.5.1                 spatstat.sparse_2.0-0        
## [179] ResidualMatrix_1.0.0          curl_4.3.1                   
## [181] ggbeeswarm_0.6.0              leiden_0.3.8                 
## [183] openssl_1.4.4                 zip_2.2.0                    
## [185] openxlsx_4.2.3                survival_3.2-11              
## [187] rmarkdown_2.8                 munsell_0.5.0                
## [189] e1071_1.7-7                   GenomeInfoDbData_1.2.4       
## [191] haven_2.4.1                   reshape2_1.4.4               
## [193] gtable_0.3.0                  spatstat.core_2.1-2
import session_info
session_info.show()
## -----
## anndata             0.7.4
## numpy               1.19.5
## scanpy              1.6.0
## scipy               1.5.2
## session_info        1.0.0
## -----
## Python 3.7.9 (default, Aug 31 2020, 07:24:53) [Clang 10.0.0 ]
## Darwin-18.7.0-x86_64-i386-64bit
## -----
## Session information updated at 2021-08-02 18:47